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ABSTRACT 

Multi-dimensional upwind-differencing schemes for the Euler equations are reviewed. On 
the basis of the first-order upwind scheme for a one-dimensional convection equation the two 
approaches to upwind differencing are discussed: the fluctuation approach and the finite- 
volume approach. The usual extension of the finite- volume method to the multi-dimensional 
Euler equations is not entirely satisfactory, because the direction of wave propagation is 
always assumed to be normal to the cell faces. This leads to smearing of shock and shear 
waves when these are not grid-aligned. Multi-directional methods, in which upwind-biased 
fluxes are computed in a frame aligned with a dominant wave, overcome this problem, but 
at the expense of robustness. The same is true for the schemes incorporating a multi- 
dimensional wave model not based on multi dimensional data but on an ’’educated guess” 
of what they could be. 

The fluctuation approach offers the best possibilities for the development of genuinely 
multi-dimensional upwind schemes. Three building blocks are needed for such schemes: 
a wave model, a way to achieve conservation, and a compact convection scheme. Recent 
advances in each of these components are discussed; putting them all together is the present 
focus of a worldwide research effort. Some numerical results are presented, illustrating the 
potential of the new multi-dimensional schemes. 
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1 Introduction 


GFD algorithms for the coming generation of massively parallel computers will have 
to be extremely robust. They will most likely be implemented on adaptive unstruc- 
tured grids, and will be used for ambitious simulations of steady and unsteady three- 
dimensional flows. In such a complex environment there is little place left for hand- 
tuning parameters that regulate accuracy, stability and convergence of the computa- 
tions. A typical algorithm will make very intensive use of local data, with a minimum 
of message passing. 

Algorithms of this nature exist already in CFD: they are the upwind-differencing 
schemes, computationally intensive but unsurpassed in their combination of accuracy 
and robustness. While these favorable properties are explainable for one-dimensional 
methods, it is a stroke of luck that upwind schemes work as well as they do for 
two- and three-dimensional flow. Their design is commonly based on one-dimensional 
physics, namely, the solution of the one-dimensional Riemann problem that describes 
the interaction of two fluid cells by finite-amplitude waves moving normal to their 
interface. The inadequacy of this technique clearly shows up when the numerical 
solution contains shock or shear waves not aligned with the grid, for instance, by a 
loss of resolution. 

The need to incorporate genuinely multi-dimensional physics in upwind algorithms 
was recognized as early as 1983 by Phil Roe [1], A study of discrete multi-dimensional 
wave models by Roe followed in 1985 (ICASE Report 85-18, also [2]), but it took until 
1991 [3] before any algorithms based on such wave models became truly successful. 
Important contributions to this development were made by Herman Deconinck and 
collaborators [3, 4] at the Von Karman Institute in Brussels. The new upwind schemes 
are formulated on unstructured grids with data in the vertices of triangular or tetra- 
hedral cells. 

While genuinely multi-dimensional methods were slowly developing, partial suc- 
cesses were booked by putting some multi-dimensional information into the Riemann 
solvers used in conventional upwind schemes. In particular, it became the fashion 
to obtain a plausible wave-propagation angle from the data, rather than accepting 
the angle dictated by the grid geometry. The earliest work of this kind is due to 
Steve Davis [5]; it recently was picked up by a number of authors: Levy, Powell and 
Van Leer [6], [7], Dadone and Grossman [8, 9], Obayashi and Goorjian [10], Tamura 
and Fujii [11]. Roughly speaking, they apply Riemann solvers in several, physically 
appealing, directions; I shall refer to their work as the multi-directional approach. 

Related, but closer to the genuinely multi-dimensional approach is the work of 
Rumsey, Van Leer and Roe [12, 13, 14, 15] and Parpiaand Michalek [16, 17]. These au- 
thors independently developed almost identical multi-dimensional wave models based 
on minimizing wave strengths. These wave models requires only two input states, just 
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Figure 1: Two views on scalar upwind differencing: (a) nodal-point interpretation; 
(b) finite- volume interpretation. 

as a regular Riemann solver. 

In support of these quasi-multi-dimensional approaches, aimed at putting better 
physics into interface fluxes, some authors have dedicated efforts to improving the 
interpolation or reconstruction step that precedes the flux calculation. On a struc- 
tured grid the reconstruction of a non-oscillatory distribution of flow variables from 
their cell-averages usually is done dimension by dimension; a fully multi-dimensional 
reconstrucion is indispensable in achieving higher accuracy. Barth and Frederickson 
[18] indicated how to reconstruct a smooth function up to arbitrarily high order on 
an unstructured triangulation; Abgrall [19] showed how to implement truly multi- 
dimensional limiting of higher derivatives. 

In this lecture I shall review a decade of efforts toward multi-dimensional upwind- 
differencing, with the accent on the very latest developments. The discussion is limited 
to the multi-dimensional physics that goes into these methods; multi-dimensional 
reconstruction will not be further mentioned. For a somewhat different emphasis 
or point of view the reader is referred to three excellent other reviews of multi- 
dimensional methods [20, 21, 4] that have been presented in the past year. 


2 Two views of one- dimensional upwinding 


In order to appreciate the problems surrounding multi-dimensional upwinding it is 
useful to consider the principles of one-dimensional upwinding. The reader is assumed 
to be familiar with the theory of conservative upwind schemes; as a tutorial Roe’s [22] 
review article is recommended. 

Upwind differencing is a way of differencing convection terms. For the scalar 
convection equation 

u t + cu x = 0, (1) 

the simplest upwind-difference scheme, of first-order accuracy, reads 
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Figure 2: Two approaches to upwind differencing for the Euler equations: (a) fluctu- 
ation approach; (b) finite-volume approach. 
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Scheme (2,3) can be regarded as a formula for updating, from t n to £ n+1 , either the 
nodal-point value of u in or the cell average of u in cell i. These two view-points 
are illustrated in Figures la and lb. The distinction is significant, because it leads 
to distinct methods for more complex equations. In the development of schemes for 
the one-dimensional Euler equations, the first view-point has led to the concept of 
fluctuation splitting , due to Roe [23, 22]; the second view-point is that of Godunov 
[24] and has led to the projection/evolution or reconstruction/evolution concept of 
finite-volume schemes, due to Van Leer [25, 26, 27]. Below I shall review the formulas 
pertinent to each approach. 


2.1 Fluctuation splitting 

Assume the system 

U t + F{U) X = 0 (4) 

represents the Euler equations in conservation form, i.e., U = (p, pu, pE) T is the 
vector of conserved state quantities and F(U) = {pu, pu 2 +p,puH) T is the vector of 
their fluxes. The equation shows that any local imbalance of the fluxes causes the 
local solution to change in time. Such a local imbalance is called a fluctuation by Roe 
[28, 1]. If source terms are present, their value must be included in the fluctuation 

[22 j. 

Define the matrix A(U) as the derivative of F(U) with respect to U, so that 

dF{U) = A(U)dU. (5) 

It is essential for the technique of fluctuation splitting that this differential relation 
be replaced by an exact finite-difference analogue, namely, 

A F = AA(7, (6) 

where A indicates a difference between neighboring nodal points. Roe [29] has in- 
dicated how to construct a mean value A of A such that Eq. (6) holds exactly for 
arbitrary pairs of state vectors. For a calorically perfect gas a suitable mean value 
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can easily be obtained by introducing the the parameter vector to = yfp ( 1 , u, H) T . 
Since both and F(w) are quadratic in the components of w : it follows that Eq. 

(6) is satisfied by A — A (fy(rD)), where w is the algebraic average of w. 

Fluctuation splitting requires that the matrix A be split into its positive and 
negative parts, i.e., 

A = A + + A-. (7) 

so that 

A F = A + AU + A~AU. (8) 

A popular name for this procedure is ”flux-difference splitting”; the term “fluctuation 
splitting is preferable because it includes source-term splitting. The first term on the 
right-hand side combines disturbances that propagate forward; in consequence, this 
tei m is used to update the right nodal point. The second term combines backward- 
moving disturbances and is used to update the left nodal point. This concept is 
illustrated in Figure 2a. Conservation is ensured because the two terms add up to a 
perfect flux difference. The first-order update formula becomes 
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In practice it often pays to abandon the matrix notation and expand A U and A F 
in terms of the individual disturbances. This yields 


3 


A U = 

Aj=1 

3 

(10) 

A F = yA k a k R k , 

(11) 


where At is an eigenvalue of A , is the corresponding eigenvector, and a*, is the 

wave strength; note that Eqs. (6) and (10) imply Eq. (11). By considering that each 
fluctuation may move forward or backward through the grid, we recover the splitting 
formula (8): 


= y a k&kRk + y a k&kRk 

AfcCO A fc >0 
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= y + y \pcY k R k (12) 

A-l A-=l 

= A + AU + A~AU. 

2.2 Finite-volume approach 

In the finite-volume approach the focus is on the numerical flux function F (Ui, Ur), 
a recipe for computing the interface fluxes from the states Ul and Ur on the left 
and right sides of the interface. The generic formula for updating cell averages of the 
conserved quantities is 


up +l = u n 


Ax 



(13) 


4 



In Godunov's first-orcler scheme the interface flux is taken from the solution at 
t > 0 of Riemann’s initial- value problem with input data 

U{x,0) = U L , x > 0, 

U(x, 0) = Ur, x < 0; 

this is illustrated in Figure ‘2b. 

For many applications it is not necessary to use the exact solution to this problem, 
hence the activity in the research area of “approximate Riemann solvers [23, 30, 31]. 
Adopting Roe’s [23] approximate solution, which is the exact solution of the locally 
linearized equation 

U t + AU x = 0, (16) 

we find three equivalent formulas for the interface flux: 

F{U l ,Ur) = F l + A~AU, 

F{U l ,Ur) = Fr-A + AU, 

F{U l ,Ur) = { -{F l + F r )-\A\AU, 

where 

\A\ = A + - A~. (20) 

In practice the formula (19) is preferred because of its symmetry; the expanded form 
is 3 

F(U l , Ur) = Uf l + Fr)- 1 -J2 \h\<x k Rk. (21) 

- Z k=\ 

Inserting the flux (19) into the finite- volume scheme (13) yields an scheme that, 
with the help of the identity (6), reduces precisely to the fluctuation scheme (9). Yet, 
there exists an important difference between Eqs. (19,13) and Eq. (9): in the latter 
the matrix A must satisfy the identity (6) in order to maintain conservation, while 
in Eq. (19) the matrix \A\ may be derived from any average A without endangering 
conservation. The flux formula (19), due to Van Leer [32, 33], preceded the fluctuation 
approach of Roe [23], based on (6), by a decade. 
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3 Intermezzo: how good is one-dimensional up- 
winding? 

To appreciate the superior accuracy and robustness of upwind differencing in one 
dimension, consider the numerical results shown in Figure 3 and 4, taken from [34] 
and [35], respectively. In Figure 3a the exact and discrete Mach-number distributions 
for choked flow through a converging-diverging channel are superimposed. First- 
order fluctuation splitting was used, including source-term splitting [22, 36] and a 
special splitting near the sonic point [34]. Although the update formula is only first- 
order accurate, it can be shown that the scheme yields second-order accurate steady 
solutions. In fact, in the steady state the scheme reduces to the two-point box scheme 
on all meshes except near a sonic point and inside a shock structure, where it becomes 
a three-point scheme. This yields the smooth transition through the sonic point and 



the crisp shock transition in the displayed results. Figure 3b shows the residual- 
convergence histories for three increasingly powerful! marching techniques: global 
time-stepping, local time-stepping and characteristic time-stepping [35]; these look 
uneventful. In Figure 4a a shockless transonic solution is reached from initial values 
containing 7 shocks and 8 sonic points; again, the residual-convergence history in 
Figure 4b for local time-stepping shows nothing unusual. 

It is this type of performance we wish to preserve when extending upwind differ- 
encing to higher dimensions. 



Figure 3: Choked flow through a converging-diverging channel, computed with a 
fluctuation scheme, (a) Initial and final Mach-number distributions; (b) residual- 
convergence histories for global, local and characteristic time-stepping. 



Figure 4: Transonic flow in a converging-diverging channel, computed with a fluctua- 
tion scheme, (a) Initial and final Mach-number distributions-; (b) residual-convergence 
history for local time-stepping. 
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Figure 5: Extending the finite- volume method to two dimensions by solving one- 
dimensional Riemann problems at all cell faces. The arrows symbolize the exchange 
of information between cells in the direction normal to their interface. 

4 Multi-dimensional extension of the finite- volume 
method 

The standard way to extend upwind differencing to the multi-dimensional Euler equa- 
tions is still the same as indicated by Godunov et al. [37] in 1961. For first-order 
accuracy, initial values are assumed to be constant in each cell, just as in one di- 
mension; fluxes at cell interfaces again follow from solving one-dimensional Riemann 
problems of the type (14,15), with x now measuring distance along the normal to the 
interface. This is illustrated by Figure 5. 

It is the projection of the true initial values onto cellwise constant distributions 
(or linear [25, 26] or quadratic [25, 38, 39] or even higher-order distributions [40]) that 
creates discontinuities at the interfaces. This leads us to introducing plane wave fronts 
parallel to the interface, and selecting, out of all possible directions, the interface 
normal as the direction for wave propagation. If the solution contains only shock 
and/or shear waves aligned or nearly aligned with the grid, this choice happens to be 
the correct one, and high resolution of such waves can be achieved in the steady state, 
just as in one dimension. Ifi however, such waves are far from aligned with the grid, 
they get misrepresented by the upwind scheme as pairs of grid-aligned waves, as shown 
in Figure 6 for a shear wave. Thus, a grid-oblique stationary wave may be represented 
by several grid-aligned running waves, leading to higher numerical dissipation and a 
considerable loss of resolution. 

Another purely numerical artifact caused by grid-aligned upwinding is the pres- 
ence of pressure disturbances across a grid-oblique shear layer. First observed by 
Venkatakrishnan [41], the explanation was provided by Rumsey et al. [12]; this phe- 
nomenon is further discussed in Section 4.2. 

From the above critique one should not conclude that in higher dimensions the 
standard upw r ind methods are inferior to other methods; the loss of accuracy just is 
much more obvious for upwind methods. 
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compression shear 


Figure 6: Misinterpretation of a grid-oblique shear wave by grid-aligned upwinding. 



Figure 7: Fluxes in a frame aligned with a wave front oblique to the grid lines. 

4.1 Multi-directional methods 

The smearing of oblique shock waves in numerical solutions has received considerable 
attention, and a proportionally large research effort has been spent in mending this 
weakness. The prevailing idea is to solve the Riemann problem in a direction more 
appropriate than the grid direction. One immediate consequence of leaving the grid- 
aligned frame is that solving one Riemann problem no longer suffices. Figure 7 shows 
that, in two dimensions, both flux vectors in the rotated frame are needed for the 
construction of the fluxes normal to the interface. 

Consider, for example, Figure 8, showing a rotated coordinate system aligned with 
level lines representing a shock front in a discrete solution. It makes sense to solve 
a one-dimensional Riemann problem in the direction normal to the front, i.e. using 
the flow-velocity components in that direction; this yields the flux in the normal 
direction. The input states for the Riemann solver are Ul± = Uj_, and { J R1 = Ur. The 
flux tangential to the shock should be obtained from state values located at Ly and 
R\\- using Ui and Ur once more would completely destroy the effect of the rotation 
[7, 14]). These values could be approximated by 

Ul\\ = Ur\\ — ^ {Ul + Ur) i (22) 

this, however, implies central differencing along the shock and leads to odd-even 
decoupling in that direction [6, 7, 42]. 
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Figure 8: A simple multi-directional flux formula. 



Figure 9: Input states for the Riemann problems in the flux computation according 
to Levy et al. 

In the work of Davis [5, 43], dating back as far as 1983, the computation of the 
tangent flux actually is more complicated than that of the normal flux. The more 
recent work of Levy et al. [6, 7, 42] and Dadone and Grossman [8, 9] is more 
mature in that the fluxes are^reated without distinction. Figure 9 shows how pairs 
of input states to the two Riemann problems, (Ull, Urx) and (U L y, Ur\\), are selected 
according to Levy et al. In their first-order method the input states in the rotated 
frame are obtained by linear interpolation between neighboring states in a ring of cells 
surrounding the interface; Dadone and Grossman simply take the value in the nearest 
cell, which apparently adds to the robustness of the method. Another, wider ring of 
cells is needed for achieving second-order accuracy. 

Various choices can be made for the rotation angle of the frame in which the 
Riemann problems are solved. A sensitive quantity is the direction of the velocity- 
difference vector, Vr — Vl-i which was adopted by Davis and also is crucial to the 
approach of Rumsey and Parpia (see Section 4.2). Levy et al. use the direction of 
the velocity-magnitude gradient V|V|, which can detect both shock and shear waves, 
while Dadone and Grossman use the pressure gradient Vp, which only detects shocks. 
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For a more detailed description of the nnilti-directional approach the reader mav 
be referred to reference [9] in these proceedings. 

After a decade of multi-directional methods, what benefits have been demon- 
strated? Surely, these methods yield impressive results when applied to first-order 
schemes: shock and shear waves not aligned with the grid are represented as if 
computed with a higher-order method. The improvement brought to higher-order 
schemes, though, is a lot less spectacular, and this is understandable. On the one 
hand, there is not much room left for a further reduction of wave spread (more for 
shear waves than for shock waves); on the other hand, loss of monotonicity may occur, 
against which there are no effective limiters, and convergence to a steady state suffers 
under the strong nonlinearity of the methods. 

In my opinion, the multi-directional approach has had a clear impact on computa- 
tional fluid dynamics. Although complete multi-directional methods will survive only 
if the problem of ensuring robustness can be solved, I expect that elements of such 
methods may find their way into standard, direction-split codes, to help resolve flow 
features arising in specific flow problems. 

4,2 Minimum-strength wave models 

In the work of Rumsey, Roe and Van Leer [14] and Parpia and Michalek [17], the 
orientation of the cell interface is de-emphasized. The spatial discretization is no 
longer regarded as generating a discontinuity along the interfaces; instead, an attempt 
is made to find out what waves are actually propagating near the interface. This, of 
course, requires data spanning a multi-dimensional part of space; if only the two 
states U i and Ur are to be used, a theoretical conjecture must make up for the 
missing information. 

In the basic wave model of Rumsey et al. a special set of 4 waves is used to 
match the state difference Ur — Ui ; for uniqueness, the sum of the wave strengths 
is minimized. Three of these waves follow from solving a one-dimensional Riemann 
problem in the direction of the velocity difference AV, the fourth wave is a shear 
w r ave normal to the other three. This choice of waves makes sense from a kinematic 
point of view, as illustrated by Figure 10. It shows that a velocity difference AV can 
be explained by an acoustic wave traveling in the direction of AF as well as a shear 
wave traveling in the normal direction. Which explanation is the more likely one 
may be determined by also considering the pressure difference ])r — p a large value 
favors the acoustic explanation, while a small value favors the shear explanation. The 
minimization procedure takes the full state difference Ur — Ui and comes up with a 
plausible explanation in terms of all four waves. The method of Parpia and Michalek 
differs only in the choice of the functional that is minimized. Figure 11 shows the 
configuration of the plane waves crossing the interface. In practice both methods 
include a fifth wave, a weak shear wave, which corrects for the difference between the 
true direction of AV and the direction actually used; the latter may have been held 
over from a previous iteration (’’frozen”), for improvement of convergence. 

The word ’’plausible” used above indicates that the minimization procedure only 
makes an educated guess: it is possible to compose a set of initial values that is totally 
misinterpreted. Consider, for instance, the head-on collision of two gases that have 
equal, negligible pressures. In reality two strong shocks are formed, moving into the 
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Figure 10: Shock or shear? 

gases. The procedure sees as input a velocity difference not accompanied by a pressure 
difference, hence calls for a single shear wave, as if the gases avoided collision! 

The flux formula based on the above wave model is worth some discussion. As- 
suming the system 

Ut + F(U) x + G(U) y = 0, (23) 

with flux .Jacobians A(U) and B(U ), represents the two-dimensional Euler equations, 
we may again write A U as a sum: 

A U = £>,/?,. (24) 

k=\ 

The vector Rk is now an eigenvector of the matrix 

A cos 9k + B sin 9 k , (25) 

where 9k indicates the propagation angle of the k - th wave; the matrices A and B are 
standard Roe-averages. The upwind-biased interface flux is defined by 

1 5 

F(U L , U R ) =-(Fl + Fr) - J2 l A * COS (°k - ^normal) | O k R k , ( 26 ) 

2 k= 1 

i.e. still by formula (21), but with the wave speeds A*, projected onto the interface 
normal. Although this formula seems trivial, it must be pointed out that there no 
longer exists a relation between A F and Af/ like (6). 

In numerical practice minimum-strength wave models appear to bring the same 
benefits and problems as multi-directional methods: great improvements in shock and 
shear resolution for first-order methods, much smaller improvements for second-order 
methods, and possible loss of monotonicity and convergence. 

To illustrate the performance of this class of methods, consider Figures 12a and 
l‘2b. Both show pressure plots for steady viscous flow over a NACA 0012 airfoil at 
3° angle of attack and Reynolds number 5000, computed on a 129 x 49 O-grid by 
Rumsey [12, 14]. Under these conditions the flow separates from the upper surface, 
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Figure 11: Plane waves crossing a cell face according to the model of Rumsey et al. 

producing a detached shear layer oblique to the grid. For the results of Figure l‘2a a 
second-order MUSCL-type scheme [26, 44] was used, with Roe’s [23] standard grid- 
aligned Riemann solver. The Riemann solver misinterprets the oblique shear as an 
grid-aligned shear plus an acoustic wave (see Figure 6); the latter causes a pressure 
rise or drop at the interface. Correspondingly, the steady solution shows pressure 
fluctuations across the shear layer, so that its presence can actually be detected in 
pressure plots. A grid-refinement study shows that the disturbances scale with the 
mesh size. This phenomenon was first observed by Venkatakrishnan [41] and correctly 
explained by Roe; in fact, it motivated the work of Rumsey, Van Leer and Roe. As 
seen from Figure 12b, the minimum-strength wave model properly recognizes the 
oblique shear layer and generates clean pressure contours. 

The same method gives an unexpected improvement in the representation of in- 
viscid stagnating flow. The explanation is found in Figure 13, showing the turning 
of the flow near a stagnation point S as represented by the discrete velocities in the 
three cells marked 1, 2 and 3.. A grid-aligned Riemann solver interprets the velocity 
difference between vertical neighbors 1 and 2 as a compression (V^i > V y 2 ), and the 
velocity difference between horizontal neighbors 2 and 3 as an expansion (l 4-2 < K*); 
this leads to pressure variations of the order of AVA The wave model detects only 
very small pressure changes (A p ~ pA(V 2 )) and therefore explains both velocity dif- 
ferences by shear waves. Although this still is not the right explanation, the result is 
a decrease in numerical entropy production. The effect is rather large for first-order 
methods, as can be judged from Figure 14 showing entropy contours for inviscid flow 
over a NACA 0012 airfoil at M — 0.3, a = 1°, on a sequence of O-grids. The reduced 
entropy levels lead directly to reduced numerical drag levels, as Figure 15a shows. For 
second-order schemes the effect, as usual, is less dramatic; the drag values are given 
in Figure 15b. 
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Figure 12: Viscous separating flow over a NACA 0012 airfoil at M = 0.5, a — 3° and 
Re = 5000. Pressure contours on a 129 x 49 C-grid, obtained with a second-order 
upwind scheme incorporating (a) Roe’s grid-aligned Riemann solver; (b) the five-wave 
model of Rumsey et al. 


5 Multi-dimensional fluctuation approach 

The fluctuation approach to upwind differencing lends itself better to extension into 
higher dimensions than the finite-volume approach. Recall that a fluctuation is a 
local flux imbalance causing a non-zero time derivative of the local solution. For the 
one-dimensional Euler equations (4) the quantity -A F equals the residual evaluated 
on a one-dimensional mesh: 


/ U t dx = - [ F x dx = -A F. 

J mesh J mesh 


(27) 


This suggests extension of the fluctuation approach beyond one dimension by regard- 
ing each multi-dimensional mesh residual as the sum of a finite number of waves (say, 
rn), moving in all possible directions. Thus we discretize the two-dimensional Euler 
equations as 


/ / (Jtdxdy = - <f(Fdy - Gdx) = ^ \ k a k R k , 

J ./mesh J k~\ 


(28) 
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Figure 13: Turning of the flow in three cells near a stagnation point 5 at a wall. 




129 x 37 




(a) (b) 


Figure 14: Entropy contours for inviscid flow over a NACA 0012 airfoil at M = 0.3, 
a - 1°, generated on a sequence of O-grids with a first-order scheme incorporating 
(a) Roe’s grid-aligned Riemann solver; (b) the five- wave model of Rumsey et al. 

where the matrices A and B are multi-dimensional averages that remain to be de- 
fined. Since the fluctuation approach is a nodal-point approach, and we wish to 
develop only schemes of maximum compactness, we shall use a grid of triangular 
meshes, with data given in the nodal points. For the computation of the residual 
on such meshes it suffices to apply the trapezoidal integration rule on each side of 
the triangle. The fluctuations resulting from residual decomposition must be sent to 
the triangle’s vertices according to some distribution scheme that approximates the 
convection equation. 

It follows that, for the construction of a genuinely multi-dimensional upwmd- 
differencing scheme, three components are needed: 

1 . A reliable multi-dimensional wave model for representing the residual; 

2. A way to ensure conservation, i.e. a multi-dimensional extension of Roe’s matrix 
average; 
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Figure 15: Grid-convergence study of the drag coefficient based on (a) the first-order 
solutions of Figure 14; (b) the corresponding second-order solutions. 

3. A multi-dimensional convection scheme for advancing the waves. 

Each of these will be discussed in a separate subsection. 


5*1 Multi-dimensional wave models 


The modeling of a local Euler residual by a finite number of waves was launched as 
a research subject by Roe [2]; his first paper, however, gave no specific instructions 
as to how the model would be used in a numerical integration of the Euler equations. 
This is not surprising, given that the other problems - multi-dimensional conservation 
and advection - had not yet been addressed. 

The latest version of Roe’s wave model calls for four acoustic waves, running along 
the principal strain axes of the local fluid element, a shear wave making a 45° angle 
with the acoustic waves, and an entropy wave running in the direction of the entropy 
gradient; see Figure 16. Thus, m = 6 in Eq. (28). These six waves are defined 
by two independent angles and six strengths; therefore, eight independent pieces of 
information need to be supplied per triangular mesh. This information is available 
in the form of the gradient of the state vector; its mesh value is computed with the 
trapezoidal rule from the following boundary integrals: 


U x = 
Uy = 



-7— l Udy\ 

AvtO, J mesh 


l 


AvtCL J mesh 


/ 

Jrn 


Uix. 


(29) 

(30) 


A detailed discussion of this wave model, including the three-dimensional case, can 
be found in Roe’s contribution to the present volume [45]; numerical results obtained 
with this model are presented in the contribution by Catalano et al. [46]. 

This section would not be complete without a discussion of the work of Hirsch and 
collaborators [47, 48, 49]. Their multi-dimensional approach is based on diagonalizing 
the Euler equations, i.e. changing these into a system of convection equations, by a 
transformation of state variables. The transformation itself depends on the local 
gradient of the solution, making the diagonalization essentially nonlinear. For certain 
data the transformation does not exist, in which case it is chosen so as to minimize 
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Figure 16: Roe’s two-dimensional six-wave model. The acoustic waves run parallel 
to the principal strain axes (dashed); the strain ellips (dotted) shows the kinematic 
deformation of a circular fluid element. 

the off-diagonal terms. The update scheme, though, can be made identical to a 
fluctuation-based scheme: decomposition of the residual along certain eigenvectors, 
followed by convection of the components [50]. In two dimensions the diagonalization 
is equivalent to using one particular four-wave model; clearly, the fluctuation approach 
offers much more flexibility. 

5.2 Multi-dimensional conservation 

The multi-dimensional extension of Roe’s averaging of the flux Jacobian was indepen- 
dently discovered by Roe and Struijs, and is presented in a joint paper [51]. This very 
recent (1991) addition to the multi-dimensional toolbox applies exclusively to trian- 
gular meshes in two dimensions and tetrahedral meshes in three dimensions. The 
following description and explanation of the two-dimensional averaging apply to the 
special case of a calorically perfect gas. 

To begin with, assume that the parameter vector w = y/p{l, u, v, H) is distibuted 
linearly over a mesh triangle with vertices labeled 1, 2 and 3. Denote the average of 
w over the triangle by w\ we then have 

w = -(wi W-2 + w 3 ). (31) 

As before, U(w) and F( w), and also G(U), are quadratic in the components of w , so 
that the Jacobian matrices U w , F w and G w are linear in w, and therefore also in x 
and y. Considering that U x = U w w x , U y = U w w y , etc., where w x and w y are constant 
over the entire triangle, we conclude that Vi/, VF and VG also_vary linearly over the 
triangle. Using the definition of the mesh-averaged gradient Vi/ given in Eqs. (29), 
(30), and similar definitions of VF and VG, we easily derive the relations, 

VF = A(t/(u>))Vi/, 

^G = B(U(w))VU, 
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(32) 

(33) 




Figure 17: Stencils of two-dimensional upwind convection schemes; case a > b > 0. 

(a) Sidilkover’s second-order scheme. The fluxes for cell 1 nominally are compurted 
by linear interpolation between upstream pairs of data, but the fluxes at the North 
and South faces must be limited to prevent numerical oscillations. The limiters are 
based on the ratios a(u t - u 2 )/[b{u 5 - u 2 )] and a(u 3 - u 4 )/[b(u 2 - u 4 )], respectively. 

(b) Standard second- or third-order grid-aligned scheme. 


which are direct extensions of the one-dimensional relation (6). The extension to 
three-dimensional averaging is self-evident. 


5.3 Multi-dimensional convection 

The pursuit of multi-dimensional convection schemes has kept a number of authors 
busy over the past three years. In two dimensions the basic equation to be solved is 

u ( + au x + bu y = 0, ('^0 

where a and b are constant velocity components, or, in vector notation, 

u t + S • Vu = 0 (35) 

The first significant work was that of Sidilkover [52], who, among other things, 
showed how a second-order upwind scheme, with residual computed on a square mesh, 
can be made noil-oscillatory by standard limiters without undue spreading of the 
stencil. The domain of dependence for this algorithm is shown in Figure 17a, for the 
case a > b > 0; note how compact this is in comparison to the stencil of a standard 
second-order upwind scheme, shown in Figure 17b [27]. He also coined the name 
”N-scheme” for the first-order scheme that, on a cartesian grid, takes its data from 
the upwind triangle fitting the convection path most tightly (N stands for narrow). 
For example, for point 1 in Figure 17a it would be triangle (124). This scheme, as 
shown in [53], is optimal in the sense that, among all schemes with upwind triangular 
domain of dependence, it combines the smallest truncation error with the largest 
stable time-step. The three-dimensional extension is also described in [53]. 

While the triangles in Sidilkover’s work were still considered subdivisions of squares, 
they become autonomous in later work by other authors. A major step in the devel- 
opment of two-dimensional convection schemes was the realization that there are two 
types of triangles [54]: those with one inflow side and those with two inflow sides. 
This is illustrated in Figure 18. If there is only one inflow side, the fluctuation ap- 
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Figure 18: Two kinds of triangles: (a) with one inflow side; (b) with two inflow sides. 

proach dictates that the entire residual be used to update the opposite node. This is 
the unique ’’single-target ’ form of the scheme, similar to the one-dimensional upwind 
scheme 2. If, however, there are two inflow sides, it may be argued that the residual 
be distributed over the two nodal points defining the third side. This is the ’’dual- 
target” form of the particular scheme; each choice of distribution weights defines a 
new scheme. The spreading of the residual information over two points implies a 
potential loss of resolution, inherent to multi-dimensional numerical convection; there 
is no one-dimensional analogue of this effect. 

In the development of multi-dimensional convection schemes, three design criteria 
play a decisive role. According to these, it is desirable for a scheme to be 

1 . linear: for a given grid geometry and flow angle the solution depends linearly on 
the data. This promotes convergence to a steady numerical solution. It is well 
known that the presence of nonlinear devices in the scheme, such as limiters [44] 
and frame rotation (see Section 4.2) can slow down or even halt the convergence 
process; 

2. linearity preserving (LP): data of the form 

u(x,y) = bx - ay, (36) 

which is a steady solution of Eq. 34 are not changed by the scheme. This 
promotes the accuracy of the scheme. It can be shown [54] that LP schemes 
yield second-order-accurate steady solutions of Eq. 34; 

3. positive: the scheme has positive coefficients. This is sufficient for preventing 
numerical oscillations. 

From one-dimensional finite-difference theory we know - and have known so for 
a long time - that the above conditions are mutually exclusive. There is a famous 
theorem by Godunov [24] which says that no linear convection-diffusion scheme with 
positive coefficients can be more than first-order accurate. With reference to our de- 
sign criteria for multi-dimensional convection schemes this theorem reads: 

There are no linear positive LP schemes. 

Again, nonlinearity is essential for the design of accurate, non-oscillatory schemes. 
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Figure 19: Dual-target form of convection schemes: (a) N-scheme; (b) LDA-scheme. 

Among the various upwind convection schemes proposed in recent years, three 
schemes stand out; these are discussed below. They all are as compact as can be, 
requiring data on only one triangle for the approximation of the convection equation. 
A small miracle is that even positivity can be achieved without leaving the triangle . 
Of course, each nodal point is a vertex of a number of triangles and may receive 
fluctuations from several of these; programming therefore must be triangle-based. 
Some results of numerical experiments are presented in Section 6. 


The N-scheme: the optimal linear positive scheme 

The name of this scheme suggests equivalence to Sidilkover 's N-scheme, but it actu- 
ally is more general. Sidilkover's scheme is just the single-target form, common to 
all compact schemes; fluctuations from triangles requiring a dual-target scheme are 
ignored in the update. The dual-target form of the current N-scheme uses distribution 
weights proportional to the components of the convection speed along the two inflow 
sides, as depicted in Figure 19a. This makes the scheme optimal in the sense of having 
the largest stability range for the time-step [54]. It is also linear and positive, and 
therefore can be no more than first-order accurate. 

The NN-scheme: the optimal nonlinear positive LP scheme 

This scheme is a nonlinear variant of the N-scheme. hence the second N. The nonlinear 
procedure included in this scheme has absolutely nothing in common with the TVD- 
enforcing limiters included in one-dimensional convection schemes. It is based on 
the observation that in the convection equation (35) the component of the convection 
velocity a perpendicular to the solution gradient Vu, has no effect on u t . We therefore 
are allowed to replace a by any velocity that has the same component parallel to Vu, 
as shown in Figure 20. This component, indicated by a w , is the velocity at which 
the level lines of u normal normal to themselves, i.e. the wave speed of the local 
distribution of u . This wave speed is the smallest of all admissible convection speeds; 
it actually vanishes with the residual. We may now adopt the following strategy: it 
both a and a w call for a dual-target scheme, we replace a by a w in the N-scheme; in 
all other cases the scheme becomes or remains a single-target scheme. In the case ot 
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u=constant 


Figure 20: The NN-scheme: nonlinear single-target form. The convection speed a 
calls for a dual-target scheme but is replaced by a x , which calls for a single-target 
scheme. The wave speed a w is the component of a and ai parallel to Vu. 

Figure 20, a is replaced by a x , the nearest admissible speed yielding a single-target 
scheme. The resulting scheme does not change any nodal value if the residual vanishes, 
hence is LP, and maximizes the allowable time step. 

Numerical results indicate that the accuracy of the NN-scheme lies between first- 
and second-order; see further Section 6. 

The LDA scheme: a non-positive linear LP scheme 

This scheme is one of several low-diffusion schemes, designed for a low truncation 
error. In the dual-target form of the scheme the distribution weights are inversely 
proportional to the areas of the triangles cut from the mesh triangle by a streamline 
through the inflow vertex; see Figure 19b. This scheme is not positive, but very accu- 
rate: on a uniform grid it achieves third order accuracy, as demonstrated in Section 
6. 


The above schemes have served as the basis for convection-diffusion schemes in a 
study by Tomaich and Roe [55]. Since the diffusion operator can not be approximated 
on a single triangle, their schemes are formulated with reference to a central nodal 
point. Numerical solutions of the Smith-Hutton [56] test problem demonstrate that 
these schemes rival the best exsisting convection-diffusion schemes in accuracy. In 
addition, their way of discretizing the Laplacean is directly applicable to any of the 
disipative terms included in the Navier-Stokes equations; thus, the basis for genuinely 
multi-dimensional Navier-Stokes codes has been laid. 


6 Numerical results 

To support some of the statements made about the new, compact convection schemes 
I first show how these schemes fared in a comparative grid-refinement study by Jens 
Muller [57]. The problem is that of convection of a Gaussian distribution over a semi- 
circle; Inflow is at y = 0, x < 0, outflow at y = 0, x > 0. Four kinds of grids were 
used, of which three examples are displayed in Figure 21. Grids a and f) derive from 



The a-mesh 


The ,3-mesh 


The 7 -mesh 



Figure 21: Three grids used in the circular-convection experiments. 



Figure 22: Grid-convergence of L 2 -gttoy for convection of a Gaussian over a semicircle 
by various schemes on various grids. 

a uniform cartesian grid by adding diagonals, in grid ct those diagonals aie chosen 
that are least aligned with the convection direction, in f3 those most aligned. Giid 
7 is a irregular perturbation to / 3 , while 8 (not shown) is a minor perturbation to 
7 . Figure 22 shows the convergence of the L 2 -e rror produced by the N-, NN- and 
LDA-schemes on the different grids. From the slope of the graphs of log(error) versus 
log(mesh-width) the following conclusions can be drawn: 

1. The N-scheme is somewhat less than first-order accurate; 

2. The NN-scheme is closer to being second-order accurate than first-order accu- 
rate; 

3. The LDA-scheme is third-order accurate on a regular grid, second-order or less 
on a perturbed grid; 
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Figure 23: Contours (left) and cut along the z-axis (right) of the solution obtained 
with the LDA-scheme on a 20 x 20 /3-grid. 



Figure 24: Inviscid flow at Mach 1.4 through an inlet, computed by a genuinely 
multi-dimensional upwind Euler code. Shown are Mach-number contours and the 
unstructured grid used. 

4. All schemes decrease their error when diagonals are aligned with the flow. 

Most surprising is the achievement of third-order accuracy on regular grids, consider- 
ing the limited amount of information going into these compact schemes. Figure 23 
gives an idea of this high accuracy by showing solution contours and a cut at y = 0 
obtained with the LDA-scheme on the very coarse /3-grid of Figure 21 (Gaussian rep- 
resented on 10 meshes). Similar results obtained with three-dimensional extensions 
of the schemes can be found in Deconinck’s comprehensive review paper [4]. 

While the search for compact convection schemes continues, several authors are 
trying to put together the ingredients listed in Section 5, producing a genuinely multi- 
dimensional upwind Euler code. Advanced numerical results can be found in the 
present volume in the contribution by Catalano et al. An earlier successful calcula- 
tion of supersonic inlet flow by Struijs et al. [3] produced the Mach-contours shown 
in Figure 24; superimposed is the moderately irregular triangulation. The results 
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Figure 25: The finite-volume scheme of Barth, Powell and Parpia, based on a trian- 
gular grid and its median dual. The flux F x across one median element located in 
triangle T is computed using wave data from triangle T and fluxes from the vertices 
Li and R. For the flux F 2 across the other element in the same triangle the same 
wave data are used, but the fluxes are taken from L z and R. 

demonstrate the excellent shock-capturing ability of the NN-scheme. 


7 Finite- volume schemes revisited 

From Section 5 the reader might get the impression that genuinely multi-dimensional 
schemes can only be formulated on triangular and tetrahedral meshes, and that they 
are incompatible with the finite-volume formulation. If this were true, it would mean 
a serious restriction on the use of such schemes within the CFD community, for it 
is not at all clear that unstructured triangular or tetrahedral grids are the way of 
the future. An alternative, for instance, is offered by adaptive cartesian grids [58]. 
The emphasis on triangles in Section 5 arises from the experience that the numerical 
building blocks, e.g. Roe’s matrix averaging, take their simplest form on such meshes. 
Therefore, in developing a multi-dimensional scheme of a different format it would be 
good practice to start with the wave decomposition of residuals on an underlying 
triangular grid. 

As an example of such practice, consider the two-dimensional finite-volume Ruler 
scheme of Powell, Barth and Parpia [59], illustrated in Figure 25. The wave model 
indeed is applied to data on triangular meshes; for the update, though, a finite- volume 
scheme is chosen. Cell faces are formed from the medians of the triangles, yielding 
the so-called median-dual grid. Across each median element of the cell contour a flux 
is computed using an equation of the form (26), where L and R denote the vertices 
of the triangle side bisected by the median, and the sum includes all waves identified 
in the triangle. Note the difference with the scheme of Rumsey et ah, where the wave 
model would be based solely on Ur - U L . The resulting nonlinear scheme, applied 
to a scalar convection equation, is LP and positive and appears to be more accurate 
than the NN-scheme. 
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8 Conclusions 


The state of the art in genuinely multi-dimensional upwind differencing has made 
dramatic advances over the past three years, owing to a shift from the finite- volume 
approach to the Actuation approach. The basic ingredients for multi-dimensional 
Euler codes, i.e. wave model, conservation principle and convection scheme, are ready 
for integration, and the first numerical results look good. The coming years will yield 
many more Euler applications in two and three dimensions, further improvements in 
wave models and compact convection schemes, and extension of the approach to the 
modeling of the Navier-Stokes equations. 
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